Application of self-consistent a method to improve the performance of model 

exchange potentials 
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Self interaction error remains an impotrant problem in density functional theory. A number of 
approximations to exact exchange aimed to correct for this error while retainining computational 
efficiency had been suggested recently. We present a critical comparison between model exchange po- 
tentials generated through the application of the asymptotically-adjusted self-consistent a, AASCa, 
method and BJ effective exchange potential advanced in [A.D. Becke and E.R. Johnson, J. Chem. 
Phys. 124, 221101 (2006)] and [V.N. Staroverov, J. Chem. Phys. 129, 134103 (2008)]. In particu- 
lar we discuss their compliance with coordinate-scaling, virial and functional derivative conditions. 
We discuss the application of the AASCa method to generate the AA-BJ potential. A numerical 
comparison is carried out through the implementation of a fully-numerical diatomic molecule code 
yielding molecular virial energies and ionization potentials approximated by the energies of the 
HOMO orbitals. It is shown that some of the shortcomings of these model potentials, such as the 
non-compliance with the Levy-Perdew virial relation, may be eliminated by multiplying the response 
term by an orbital-dependent functional a, which can be simplified to a constant determined during 
the self-consistent procedure (self-consistent a). 



PACS numbers: 



I. INTRODUCTION 



An important problem arising in the application of the 
Kohn-Sham equations is that of the construction of the 
local exact exchange potential. In principle, this poten- 
tial can be exactly calculated following the procedure 
outlined in the OEP method. 1-9] However, in prac- 
tice, there are a number of difficult numerical imped- 
iments that bar the way to the realization of this ex- 
act approach [Tol - fT3l | in addition to some deeper problems 
stemming from general instabilities of Kohn-Sham poten- 
tials in finite dimensional sub-spaces. [l3| 

For this reason, particularly in recent years, much at- 
tention has been devoted to the development of alterna- 
tive local-exchange potentials which are simple to apply 
but which at the same time yield sufficiently accurate 
results. SELlH 

A few years ago, we introduced such an approach, 
which we called the "asymptotically-adjusted self- 
consistent a" (AA-SCa) method. Although the detailed 
theoretical justification for this method is presented in 
Ref. [13], we comment here on two of its characteristics. 

The first is that in the AA-SCa method the following 
model potential is postulated: 



tial related to an approximate (and arbitrary) exchange 
functional E®[p] = J d 3 rp(r)e°([p]; r) yielding the local 
exchange potential v x (r) = &E^\p]/8p; (however, more 
generally, the response term v® esp also may be modelled) . 
In Eq. ([1]), a x [{'i/ , i}] is the functional 



£ x [{V>J]-£ xLP [z; s ([p]),p] 

SxLp[w r osp ([p]),p] 



(2) 



where £J x [{-0j}] is the exact orbital expression for 
the exchange energy, and £" x lp [v ([/£>]; r), p(r)l = 
J <i 3 ri;(r)[3p(r)+r-Vp(r)]) is the Levy-Perdew expression 
for the exchan ge e nergy functional corresponding to the 
potential v(r) |25l. [26j. The local potential given by Eq. 
(P) is called "asymptotically adjusted" because while the 
first term guarantees the correct asymptotic behavior of 
— 1/r for large r the second term does not contribute in 
the asymptotic region. 

The second characteristic follows from Eq. ©: 



#xLP [«S([P]),P] 
AASCa 



a x [{^i}]-ExLP [v- 
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p(W).p] 

(3) 



v AASCa (r) = v s (r) + a x [m]v°(r) 



The variational derivative of this functional with respect 
to the Kohn-Sham orbital yields: 



(1) 



where «s( r ) is the local Slater potential and where 
rfi fv\ — _ OtOfUl-j.) i s a l oca l response poten- 



5E AASCa [{vj z }} 
Sipj 



AASCa 



(r)+Av x (r) Vito (4) 
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where A« x (r)i/>j(r) = [v xj (r) - v s (r) - a x «° esp (r)]^(r) 
is a non-local correction to the v AASCa (r) potential. It 
has been shown [24j that the contribution of Av x (r) 
to the energy is AE x [{vji}] = where A£" x [{?/>;}] = 



2 



E x [{i>i}] - BxLp[«s] - a x [{ipi}]E xLP [v° esp \. Hence, the 
omission of A« x (r) in the Kohn-Sham equation only af- 
fects the quality of the converged orbitals but it does not 
change the expression for the energy. 

The AASCa method has been applied previously to 
improve the potentials and energies of several DFT ex- 
change functionals. In particular, we have analyzed the 
improvements brought about by this method with re- 
spect to the LDA and PW91 functionals and proposed 
two models for v® esp term [24j . 

However, quite recently a very simple potential de- 
noted as the BJ effective exchange potential has been 
proposed in a heuristic way by Becke and Johnson (20|. 
This potential contains a Slater term plus a local response 
one. This work has been extended by Staroverov [23| who 
has advanced a family of model potentials which have the 
form of a Slater potential plus a response term which is 
modeled. In all these cases the energy is evaluated using 
the exact expression i? x [{-0i}] for the exchange function- 
als constructed from the N occupied orbitals which have 
self-consistently converged for the given model potential. 

In the present work, due to the similarities between 
this family of potentials and the AASCa one, we make a 
critical comparison between these potentials. We also ap- 
ply the AASCa method using the model exchange func- 
tional associated with the BJ potential to generate the 
AA-BJ one. The systems chosen for the present compar- 
ison are some selected diatomic molecules. We show that 
application of the AASCa method docs indeed bring im- 
provements, albeit slight, on both the BJ energies and 
ionization potentials. Also, we show that it yields in- 
ternuclear distances that are in excellent agreement with 
the exact Kohn-Sham x-only results. 



II. COMPARISON WITH THE BJ AND 
RELATED MODELS 

There are three formal conditions that the exact op- 
timized effective potential (OEP) for exchange must 
satisfy [13, [2^| . These are: the variational derivative con- 
dition, the virial relation, and the scaling requirement. 

Quite clearly, the AASCa local potential v A ASCa (r) is 
an approximate one and, hence, it does not satisfy all 
three of these conditions. As it is shown in Eq. (QJ, 
the variational derivative of E AASCa [{ijji}] with respect 
to yields the potential v AASCa (r) + Au x (r). Let us 
emphasize, however, that Au x (r) does not contribute 
to the exchange energy, omission of this term in the 
Kohn-Sham equation makes the local potential v AASCa 
the approximate one. With respect to the virial rela- 
tion, it follows from the definition of v AASCa (r) given 
by Eq. (jTJ) that v AASCa (r) satisfies by construction 
the Levy-Perdew virial condition. Also, it is easy to 
show using the fact that both E x [p\] = XE x [p] and 
ExLp[v(\p\]),p\] = \E xh p[v([p]),p] that a x [{V>i}], as 
defined by Eq. @, is invariant under coordinate scal- 
ing (does not depend on A). Then using E AASCa [px] = 



E*lp [vs{[px]), P\] + u*E xLP [v^gpdpx]), px] and the fact 
that SE AASCa [p x ] = XSE A ff^ a [p], and following the 
same arguments as in Ref. [28j . it can be readily shown 
that v AASCa (r) also satisfies the scaling property. 

In the case of the potentials introduced by Becke and 
Johnson (20| and Staroverov [23|, it is clear that they 
satisfy the scaling requirement. However, the fact that 
these potentials do not satisfy the Levy-Perdew condition 
implies that the virial relation in the Born-Oppenheimer 
approximation [2t| is not satisfied either (see Ref. [3(| for 
details). Moreover, as a consequence of this, "there is no 
unique choice of functional for the evaluation of the ex- 
change energy" for structure optimization or total energy 
comparison (see Ref. [3l|, where the Becke- Johnson po- 
tential was employed for band gap calculations in solids). 
Also in the case of the BJ and other model potentials 
discussed in [13, HH, the third requirement from [27l [28|. 
namely, that the model potential must correspond to the 
functional derivative of the model functional with respect 
to the density is not satisfied. 

We denote as in Ref. [23| by E conv the total energy 
computed with the exact exchange expression £ x [{^/)j}] 
using the converged orbitals and by E V i r the corre- 
sponding one which includes the Levy-Perdew expression 
E x j_,p[v x , p] for exchange. Since for the BJ and related 
potentials -E x [{i/>i}] 7^ S X Lp[u x ,p] the total energy val- 
ues E conv and E v i T are also different. The claims that 
"(£vir — -Econv) gives an indication of how close v xa is to 
the exact functional derivative" [23j or that it serves as 
an indication of the accuracy of the calculation itself (see 
Ref. [HI]) do not seem to hold generally, as there are ap- 
proximate KS exchange potentials (such as the AASCa 
one, for example), which while differing from the exact 
one, yield, nonetheless, E v - n - — E conv . 

In the present article, in order to compare the results 
obtained using the AASCa model with those coming 
from the BJ and related potentials, we provide numeri- 
cal values, in particular, for diatomic molecules. In order 
to carry out this numerical comparison, the BJ effective 
model potentials proposed in [2(| and [23[ were imple- 
mented in a fully-numerical diatomic molecule code [34| 
(due to its numerical instability stemming from "trouble- 
some ", terms, a fact that was corroborated in our test 
calculations, the gradient-corrected model [23| was not 
implemented) . 

For completeness, we include some results obtained in 
the context of the GLLB model proposed in Ref. [lij ]. 
the localized Hartree-Fock (LHF) model [HI, HH, and 
the common energy denominator approximation (CEDA) 
[HI . We compare these results with those of the B J and 
related potentials as well as with the AA-BJ ones ob- 
tained by applying the AASCa model to the BJ func- 
tional. It is shown that the AA-BJ results are only 
slightly improved due to this application. We also in- 
clude some previous AASCa model results obtained for 
the PW91 and GLLB model functionals (in par ticular of 
the AA-PW91 and AA-m2 types, see Ref. [24j|). These 
results show that all these models are closer to EXX than 
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the BJ or AA-BJ ones in the case of diatomic molecules. 



III. RESULTS AND DISCUSSION 

Table Q] shows the Hartrce-Fock total energies and 
the energy differences for KS-x- only methods calcu- 
lated when the virial relation is used (except for the 
"BJ(conv)" column). The energy differences of the ap- 
proximate methods (eight last columns) should be com- 
pared to the exact exchange (EXX) values shown in the 
second column or to the KS(x- only) values obtained by 
the iterative procedure described in Ref. (36j (shown in 
the third column), which are very close to the EXX val- 
ues. 

All approximate methods shown in Table HI except for 
the BJ potential, are exact for the singlet state of a two- 
electron system (H2 molecule). The orbital-dependent 
methods for the response potential term (GLLB, AA-m2, 
CEDA and LHF) provide very good approximations to 
the EXX energies, the CEDA and LHF values are only 2 
mHartrees higher than the corresponding EXX energies. 
The AA-PW91 potential, taken here as an example of 
an asymptotically-adjusted potential where the response 
term is modeled by a conventional PW9f DFT functional 
was found to yield a good approximation to the EXX 
energy; the largest difference is 39.4 mHartees for the F2 
molecule (as compared to the EXX difference which is 
8.6 mHartrees). 

The BJ differences still are significantly larger than 
those arising from all other approximate methods. The 
negative value in Table U shows that the corresponding 
energy is lower than the HF value, i.e. the variational 
principle is not satisfied. Large errors in the total en- 
ergy also affect significantly the calculated atomization 
energies and the predicted equilibrium geometries. 

The ionization potentials approximated by the high- 
est occupied molecular orbital (HOMO) energy are pre- 
sented in Table |TTJ The situation is similar: the GLLB, 
AA-m2, CEDA and LHF methods provide an excellent 
approximation to the EXX values. The AA-PW9f values 
are also very close to the EXX for all molecules presented 
in Table except for the N2 and FH. The BJ potential un- 
derestimates the ionization potential by an amount of 
^30-50%. The BJ potential is calculated without shift 
as suggested in 20]. By applying the shift, the HOMO 
energies will be equal to the corresponding HF values 
(but not to the OEP ones). 



Kohn-Sham potential. The asymptotically- adjusted po- 
tentials based on the DFT approximation for the re- 
sponse term (the AA-PW9f ) are invariant w.r.t. unitary 
transformation of orbitals (as is also the BJ potential); 
however, the AA-m2 and GLLB ones are not. The ad- 
vantages of the previously proposed AA-PW9f , AA-m2 
and GLLB models (and of the CEDA and LHF methods) 
are the following: (i) For these models E v i r = E conv (no- 
tations from [2J|); these models constitute examples of 
situations where although (E v - lr — E conv ) — 0, the poten- 
tial is still only approximate; (ii) the energies obtained by 
the AA, GLLB models and by the CEDA and LHF meth- 
ods satisfy the variational principle: i? HF < E OEP < 
E^P plox , which is not the case for the BJ model; (iii) the 
virial energy for the AA and for the GLLB models is an 
excellent approximation to the OEP/EXX energy. This 
is not the case for the model potential of Refs. [ItJIlsJ], in 
spite of the fact that the BJ model is in excellent agree- 
ment with the OEP exchange potential; (iv) the ioniza- 
tion potential approximated by the HOMO energy is a 
property entirely defined by the effective potential. The 
recently proposed BJ model potential fails to adequately 
describe this property, in contrast with the AA-PW9f, 
GLLB and AA-m2 models and with CEDA and LHF 
methods, where the agreement with the EXX values is, 
in most cases, excellent. 

Some of the shortcomings of the model potentials from 
Refs. dilllH may be eliminated, however, by scaling the 
model response term by the AASCa method, as it was 
done in Ref. [24| . A rigorous variational justification for 
this t ype of scaling is given in Eqs. (f 2) through (f 5) of 
Ref. [2J]. By applying this procedure to the BJ model, 
the new AA-BJ potential (for the spin-unpolarized case) 
reads 



^ AA - BJ =t>s 



a x [mwr[m]/ P , 



(•5) 



where r is the kinetic energy density. The self-consistent 
constant a x is defined by 



£x[W]-£xlpK,p] 
E x hp[y/T[{ipi}]/p, p] 



(6) 



The asymptotically-adjusted, GLLB, CEDA and LHF 
exchange potentials, as well as the BJ model potential, 
provide an excellent approximation to the exact exchange 



The energies for the new AA-BJ model (we emphasize 
that E v ir = E conv for AA-BJ) are presented in last col- 
umn of Table I. The AA-BJ energies are slightly closer 
to the EXX values than E conv values for the original BJ 
model (BJ(conv) column in Table I). Moreover, modified 
AA-BJ model eliminates the ambiguity with regard to 
the choice of functional for the evaluation of the exchange 
energy (conventional or virial) . The ionization potentials 
presented in Table II also are slightly improved in the 
AA-BJ model as compared to the original BJ values. 

In Table UHl we present some bond lengths values ob- 
tained from several different x-only methods. We do not 
include the diatomic molecule F2 as it does not bind at 
the level of an x-only approximation. The results show 
that for H2, FH and N2, the AA-BJ bond lengths coin- 
cide up to three decimals with those obtained by means 
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TABLE I: Full-numerical Hartree-Fock (HF) total energies (in a.u.) and differences between KS-x-only and HF total energies 
(in mHartrees) calculated at the experimental geometries. "Values are taken from Ref. 6 From Ref. [24|. c From Ref. [35[ . 
d From Ref. [33]. 





HF 


EXX a 


KS(x-on(„) 


AA-PW91 6 


GLLB 6 


AA-m! 1 


CEDA C 


LHF d 


BJ(vir) 


BJ(conv) 


AA-BJ 


H 2 


-1.1336 


0.0 


0.0 


0.0 


0.0 


0.0 


0.0 


0.0 


-80.6 


0.8 


0.0 


FH 


-100.0708 


2.0 


2.2 


17.0 


5.9 


3.7 






533.1 


10.4 


6.6 


N 2 


-108.9931 


5.2 


5.7 


30.4 


11.4 


8.3 


7.7 


7.3 


239.4 


9.9 


9.5 


CO 


-112.7909 


5.1 


5.6 


31.1 


12.2 


9.1 


7.6 


7.2 


323.0 


12.6 


11.5 


F 2 


-198.7722 


8.6 


9.3 


39.4 


16.8 


13.6 






1124.9 


23.1 


18.1 



TABLE II: Ionization potential approximated by the negative of the HOMO energies (in eV). "From Ref. 0. 6 Ref. [3]. "Ref- 
6}. d N 2 ^N+( 2 E 9 ). e N 2 ->N+( 2 n u ). 'From Ref. [U . 9 From Ref. 





HF 


EXX a 


AA-PW91 6 


GLLB 6 


AA-m2 i ' 


CEDA' 


LHF 9 


BJ 


AA-BJ 


H 2 


16.2 


16.2 


16.2 


16.2 


16.2 


16.2 


16.2 


10.1 


16.1 


FH 


17.7 


17.4 


13.3 


16.5 


18.2 






9.4 


11.0 


N 2 d 


17.3 


17.2 C 


12.7 


15.5 


17.0 


17.1 




9.9 


10.4 


N 2 e 


16.7 


18.1 c 


13.9 


16.2 


18.0 


18.5 




10.7 


11.2 


CO 


15.1 


14.1 


11.1 


13.7 


15.0 


15.0 


15.0 


8.4 


9.0 


F 2 


18.2 


14.5 


13.4 


15.9 


18.7 






9.3 


11.1 



of the exact Kohn-Sham yi-only method. For the case 
of CO, there is a difference of 0.002 angstroms between 
the AA-BJ and the KS(x-onfo/) result. The BJ(vir) ge- 
ometries differ from the KS(x- only) results for all four di- 
atomics presented in Table. The BJ(conv) bond lengths 
coincides with the exact KS(x- only) results for the case 
of N2, and differ for other molecules. 

IV. CONCLUSIONS 

All potentials discussed here have the same structure: 
they comprise a Slater potential plus a modeled response 
term. As a result, the computational cost for all models is 
approximately the same, except that for the AA-PW91, 
GLLB, AA-m2, CEDA, LHF cases, where an additional 



exact-exchange energy term has to be calculated. The 
same term, however, must also be calculated for the BJ 
model potential when the total energy is obtained not by 
means of the virial relation, (where it is denoted as E v { r 
in Ref. [H|), but through the exact exchange term (in 
which case, it is denoted as E conv ). The AASCa method 
is a simple procedure which permits to transform the 
BJ model potential into a new AA-BJ potential which 
satisfies the Levy-Perdew virial relation without increas- 
ing the computational cost (in fact, this procedure may 
be applied to model potentials of any structure without 
reducing computational efficiency). However, the calcu- 
lated AA-BJ values show only a slight improvement on 
the BJ ones, except for the bond lengths, which are in 
excellent agreement with the KS x-only results. 
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